#include <math.h>
#include <cstdlib>
#include <ctime>
#include <map>
#include <algorithm>
#include <Rcpp.h>
#include <vector>
#include "encode.h"

#ifdef __cplusplus
extern "C"{
  #endif
  enum{
    ROW_OK,
    ROW_NOTSUPPORTED,
    ROW_INVALIDCODE,
    ROW_INVALIDARGUMENT,
    ROW_INTERNALERROR,
    ROW_NOMEMORY
  };
  #ifdef __cplusplus
}
#endif

using namespace Rcpp;
using namespace std;

/*---------------------------------------------------------------------------for statistics--------------------------------------------------------------------------------------------------------------*/
double member=0;
unsigned int upMax=0,rightMax=0;
unsigned long code,codeSelect,neibor;
unsigned int amountSample,total;
double xLength,yLength,size,rCopy,dist;
inline double GFunCore(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, double r){
	size = ceil(r / 2);
	xLength = ceil(xr / size);
	yLength = ceil(yr / size);
    rCopy = r;
	amountSample=xfocus.size();
	total = owin.nrow();
	upMax=0,rightMax=0;
	
	multimap<unsigned long, int> geohashMap;
	multimap<unsigned long, int>::iterator pos;

	for (unsigned int i = 0; i < total; i++){
		encode(xr, xLength, size, owin(i,0),owin(i,1), &code);
		geohashMap.insert(pair<unsigned long, int>(code, i));
	}

	for (unsigned int i = 0; i < amountSample; ++i){
		r = rCopy;
		if (xfocus[i] >0 ){
			encode(xr, xLength, size, xfocus[i], yfocus[i], &codeSelect);
			geo_neighbors(&codeSelect, xLength, yLength, size, &neibor, &upMax, &rightMax, r, 0,0);

			for (unsigned int up = 0; up < upMax; ++up){
				for (unsigned int right = 0; right< rightMax; ++right){
					for (pos = geohashMap.lower_bound(neibor + right); pos != geohashMap.upper_bound(neibor + right); ++pos){
						dist = sqrt(pow(owin(pos->second, 0) - xfocus[i], 2) + pow(owin(pos->second, 1) - yfocus[i], 2));
						if (0<dist&&dist<r){
							r=dist;
						}
					}
				}
				neibor = neibor + xLength;
			}
			if (r < rCopy) {
				++member;
				xfocus[i] = -1;
			}
		}		
	}
	return(member / amountSample);
}

inline double kFunCore(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, double r){
	size = ceil(r / 2);
	xLength = ceil(xr / size);
	yLength = ceil(yr / size);
	member = 0;
	amountSample=xfocus.size();
    total = owin.nrow();

	multimap<unsigned long, int> geohashMap;
	multimap<unsigned long, int>::iterator pos;

	for (unsigned int i = 0; i < total; i++){
		encode(xr, xLength, size, owin(i,0),owin(i,1), &code);
		geohashMap.insert(pair<unsigned long, int>(code, i));
	}

	for (unsigned int i = 0; i < amountSample; ++i){
			encode(xr, xLength, size, xfocus[i], yfocus[i], &codeSelect);
			geo_neighbors(&codeSelect, xLength, yLength, size, &neibor, &upMax, &rightMax, r, 0,0);

			for (unsigned int up = 0; up < upMax; ++up){
				for (unsigned int right = 0; right< rightMax; ++right){
					for (pos = geohashMap.lower_bound(neibor + right); pos != geohashMap.upper_bound(neibor + right); ++pos){
						dist = sqrt(pow(owin(pos->second, 0) - xfocus[i], 2) + pow(owin(pos->second, 1) - yfocus[i], 2));
						if (0<dist&&dist<r){
							++member;
						}
					}
				}
				neibor = neibor + xLength;
			}
	}
	return member*xr*yr / (amountSample*total);
}

inline double gFunCore(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, double minR, double maxR, double density){
	size = ceil(maxR / 2);
	xLength = ceil(xr / size);
	yLength = ceil(yr / size);
	member = 0;total = owin.nrow();amountSample=xfocus.size();

	multimap<unsigned long, int> geohashMap;
	multimap<unsigned long, int>::iterator pos;

	for (unsigned int i = 0; i < total; i++){
		encode(xr, xLength, size,  owin(i,0),owin(i,1), &code);
		geohashMap.insert(pair<unsigned long, int>(code, i));
	}

	for (unsigned int i = 0; i < amountSample; ++i){
			encode(xr, xLength, size, xfocus[i], yfocus[i], &codeSelect);
			geo_neighbors(&codeSelect, xLength, yLength, size, &neibor, &upMax, &rightMax, maxR, 0, 0);

			for (unsigned int up = 0; up < upMax; ++up){
				for (unsigned int right = 0; right< rightMax; ++right){
					for (pos = geohashMap.lower_bound(neibor + right); pos != geohashMap.upper_bound(neibor + right); ++pos){
						dist = sqrt(pow(owin(pos->second, 0) - xfocus[i], 2) + pow(owin(pos->second, 1) - yfocus[i], 2));
						if (minR<dist&&dist<maxR){
							++member;
						}
					}
				}
				neibor = neibor + xLength;
		}
	}
	return member / (3.14159265359*(minR + maxR)*density*density);
}

//[[Rcpp::export]]
List GFun(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, NumericVector r){
	unsigned int amountR=r.size();
	vector<double> Gval;
	sort(r.begin(),r.end());
	member=0;
	for(unsigned int i=0;i<amountR;++i){
		Gval.push_back(GFunCore(owin,xr,yr,xfocus,yfocus,r[i]));
	}
	
	return List::create(
    _["radius"] = r, 
    _["Gvalues"] = Gval
  );
}

//[[Rcpp::export]]
List FFun(NumericMatrix owin, double xr, double yr, unsigned int xRange, unsigned int yRange, unsigned int amountSample, NumericVector r){
	srand(static_cast<unsigned int>(clock()));

	NumericVector xfocus(amountSample),yfocus(amountSample);
	for (unsigned int i = 0; i < amountSample; ++i){
		xfocus[i]=double(rand() % (int(xRange) + 1));
	}
	srand(static_cast<unsigned int>(clock()+100));
	for (unsigned int i = 0; i < amountSample; ++i){
		yfocus[i]=double(rand() % (int(yRange) + 1));
	}

	sort(r.begin(), r.end());
	vector<double>Fval;
	unsigned int amountR=r.size();
	member=0;
	for (unsigned int i = 0; i < amountR; ++i){
		Fval.push_back(GFunCore(owin, xr, yr, xfocus, yfocus, r[i]));
	}

	return List::create(
	_["radius"] = r, 
    _["fvalues"] = Fval
	);
}

//[[Rcpp::export]]
List kFun(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, NumericVector r){
	unsigned int amountR=r.size();
	vector<double> kval;
	for (unsigned int i = 0; i < amountR; ++i){
		kval.push_back(kFunCore(owin,xr,yr,xfocus,yfocus,r[i]));
	}
	
	return List::create(
    _["radius"] = r, 
    _["kvalues"] = kval
	);
}

//[[Rcpp::export]]
List gFun(NumericMatrix owin, double xr, double yr, NumericVector xfocus, NumericVector yfocus, NumericVector r){
	total=owin.nrow();
	unsigned int amountR=r.size();
	double density = total / (xr*yr);  double dr = 0.05*sqrt(density);
    vector<double>gval;
	
	for (unsigned int i = 0; i < amountR; ++i){
		gval.push_back(gFunCore(owin,xr,yr,xfocus,yfocus,r[i]-dr,r[i]+dr,density));
	}
	
	return List::create(
    _["radius"] = r, 
    _["kvalues"] = gval
	);
}

